Order-disorder phase change in embedded Si nano-particles 
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We investigated the relative stability of the amorphous vs crystalline nanoparticles of size ranging between 0.8 
and 1.8 nm. We found that, at variance from bulk systems, at low T small nanoparticles are amorphous and they 
undergo to an amorphous-to-crystalline phase transition at high T. On the contrary, large nanoparticles recover 
the bulk-like behavior: crystalline at low T and amorphous at high T. We also investigated the structure of 
crystalline nanoparticles, providing evidence that they are formed by an ordered core surrounded by a disordered 
periphery. Furthermore, we also provide evidence that the details of the structure of the crystalline core depend 
on the size of the nanoparticle 

PACS numbers: 64.70.Nd, 61.46.Hk 



I. INTRODUCTION 

Nano- scale systems behave differently than ordinary bulk 
materials since, among other reasons, their physico-chemical 
properties do depend upon their size and shape. Considerable 
effort is ongoing to understand, design, fabricate, and manipu- 
late materials at such a short length scale, so as to get tailored 
properties. In particular, the identification of how the struc- 
tural features depend upon the actual thermodynamic condi- 
tions is attracting an increasing interest as it paves the way to- 
ward explaining the structure-property relationship, an issue 
of large technological impact. Among the nano-sized systems 
of current interest, silicon nano-particles embedded in amor- 
phous SiC>2 are especially important for their possible applica- 
tion as photo-emitting materials for optoelectronics, 1 2 as well 
as for the light harvesting component of solar cells™ 

A feature strongly affecting the properties of nano-sized 
semiconductor particles is whether they are crystalline or 
amorphous. In particular, it has been experimentally observed 
that the photoluminescence of Si nano-particles embedded 
in silica strongly depends (both in wavelength and intensity) 
on their crystallinity. Their structural evolution has been ac- 
cordingly characterized: Si nano-particles are initially formed 
amorphous and then transformed into crystalline upon ther- 
mal annealing at high temperatures (typically at 1100°C or 
above) PKSl During annealing, another phenomenon has been 
nevertheless observed, namely: the growth of nano-particles, 
which makes it difficult to unambiguously identify the actual 
atomistic mechanisms driving the observed evolution. The 
aim of this paper is to resolve this ambiguity by elaborating 
a thorough atomistic explanation of the observed microstruc- 
ture evolution of an embedded Si nano-particle, through com- 
puter experiments addressed to measuring its free energy in 
different states of aggregation. The main output of the present 
investigation is that we identify the most stable phase of a Si 
quantum dot as a function of its size and the thermal condi- 
tions. This result provides evidence that at the nano-scale the 
relative stability of the ordered and disordered phases could 
be otherwise than in the bulk samples. In particular, we show 
that this result is able to explain the experimental findings on 



the mechanism of formation of crystalline nano-particles jHESI 
We also investigated the atomistic structure of the states cor- 
responding to the minima of free energy, discovering that the 
ordered states are not simple crystal-like clusters; rather, they 
are made by a crystalline core surrounded by a disordered 
shell. We also found that the details of the atomistic structure 
of the crystalline core depend on the size of the nanoparticle. 
The analysis performed in this investigation is very general 
since it is addressed to any metastable state identified by atom- 
istic simulations and to characterize the order-disorder phase 
change as a function of the size of the nanoparticle and the 
temperature of the system. 

The article is organized as follows: in Sec. [IT] we describe 
the computational methods we used for the present free energy 
calculation and describe the simulation setup. In Sec. [Ill] we 



present and discuss our results. Finally, in Sec. IV we draw 
some conclusions. 



II. THEORETICAL BACKGROUND, COMPUTATIONAL 
METHODS AND SIMULATION SETUP 



This section is divided into three Subsections. In Sec. ( |II A| ) 
we introduce the technique used to calculate the free energy 
of a computational sample in a given state. This technique re- 
quires the introduction of suitable collective variables describ- 
ing the state of the system, which are described in full detail 



in Sec. ( |IIB[ ). Finally, in Sec. (HQ we present the simulation 
setup. 



A. Free energy calculation 

In this work we compute the free energy by numerically in- 
tegrating the gradient of the free energy, which is evaluated ac- 
cording to the restraint method introduced by Maragliano and 
Vanden-Eijnden. 11 This approach allows to efficiently com- 
pute the (relative) free energy in a system containing multiple 
metastable states separated by significant free energy barriers. 

Let (x) be a suitable order parameter describing the state 
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of the system, where x is the 3N vector of the atomic positions. 
Consider the following Hamiltonian: 

H k (p,x) = K(p)+V(x) + k -{d{x) - 9*) 2 (1) 

where p is the 3N vector of the atomic momenta, K(p) is the 
kinetic energy, and V(x) is the interatomic potential. 0* is a 
possible realization of the collective variable (x) and k is a 
tunable parameter. Below, we show that in the limit k —> oo 
the observable < k(6(x) - 0*) > Hk , where < ••• > Hjc indi- 
cates the average over the canonical ensemble associated to 
the Hamiltonian /4(p,x), is the derivative of the free energy 
with respect to the parameter identifying the state of the sys- 
tem at the value 6 = 6* (dF(6)/d6\ e =e* )• By simple algebra 
it can be shown that 

<k(e(x)-e*)> Hk = (2) 

d jS ~ 1 log / <ixJpexp[-j3/4(x,p)] 

~d¥ iT 

where 3f = / dxdpexp[— j3//(x,p)] is the canonical partition 
function of the system and j3 =ksT (ks is the Boltzmann con- 
stant). Since 3? is -independent its introduction does not 
affect our argument but it is necessary for the probabilistic 
interpretation of < k(6(x) — 0*) >H k > Let us now consider 
/ dxdpexp[— j3/4(x, p)]/^ in the limit of large k: 

fdxdyexp[-PH k (x,p)} 

£s x = (3) 

/ dxdpexp[-/3#(x,p)]S(0(x) - 0*) 

& 

The r.h.s of Eq. [3] is, by definition, the probability density 
Pq(0*) to find the system in a state corresponding to 0(x) = 
0*. The relation between the free energy F(6*) and the above 
probability density function is F(0*) = -j5~ l logP e (6*). As 
a result, in the limit mentioned above, Eq.[2]reads: 

lim <k(0(x)-0*) >H k =dF(6)/dO\e=e* (4) 

k— 

By (numerically) integrating the so computed 
dF(6)/d6\o=o* we straightforwardly get the F(6) curve. 

In practice, we can compute an approximation to the deriva- 
tive of the free energy on the collective variable by cal- 
culating < k(G(x) - 0*) > Hk for large k by Molecular Dy- 
namics (MD). In this case, we replace the ensemble aver- 
age < k(6(x) — 0*) >H k by a time average of the operator 
k(6(x(t)) — 6*) along the trajectory of a constant temperature 
MD in which the atomic forces are obtained from the poten- 
tial U(x) = V(x)+k/2(0(x) - 0*) 2 . Other methods could 
also be used to compute the free energy (e.g. the Blue Moon 
ensemble 12 or the umbrella sampling 13 ). The advantage of the 
method used in this paper with respect to the Blue Moon sam- 
pling is that it does not require any un-biasing, which might be 



difficult to perform, depending on the selected collective vari- 
able; on the other hand, with respect to the umbrella sampling 
the advantage is that it does not require a technique such as the 
Weighted Histogram Analysis Method (WHAM) 14 to recon- 
struct the free energy curve; rather, this quantity is obtained 
by performing a simpler numerical integration. Nevertheless, 
the accurate calculation of the free energy F(6) might require 
the calculation of < k(0(x) — 0*) >n k in more 0* points in 
comparison to the corresponding number of positions of the 
umbrella potential (i.e. number of umbrella sampling runs). 
It is worth to remark that the present method allows to com- 
pute a free energy including both the configurational and vi- 
brational entropic contributions with no approximation (apart 
those connected with finite time simulation). 

In the present investigation we need to extend the ap- 
proach described above to the case of two collective vari- 
ables, one controlling the degree of order of the nanoparti- 
cle (and, therefore, the amorphous vs crystalline phase), the 
other monitoring its size. They are both described in detail 
in Sec. MB I The use of two collective variables is motivated 
by the need of computing the relative free energy of the or- 
dered and disordered phase at a given size of the nanopar- 
ticle. The probability density function associated to such a 
free energy is the conditional probability density function to 
observe = 0* (the first collective variable) given = 0* 
(the second one), hereafter indicated as P(0*|0*). This quan- 
tity can be computed by numerically integrating the observ- 
able < k(6(x) — 0*) >H kk n where the Hamiltonian H^y = 

K(p)+V(x) + §(0(x) - 0*) 2 + f (0(x) - 0*) 2 . As explained 
before, this observable can be computed by MD. 

The procedure described above assumes that, apart from the 
collective variable 0(x) and 0(x), the system is ergodic. How- 
ever, there might be cases in which other slow variables are 
present in the system and, therefore, the calculation of the ob- 
servable < k(0(x) — 0*) >H kk f cannot be efficiently computed 
by a straightforward MD simulation. In Sec. [ill] we show 
that this is indeed the case in the present investigation. In or- 
der to overcome this problem we combine the restrained MD 
method described above with the Replica Exchange method.EI 
A similar approach has already been used in the simulation 
of rare events in which the replica exchange method has 
been used in combination with the umbrella sampling or 
Metadynamics 17 The Replica Exchange technique consists in 
running several MD simulations at different temperatures in 
parallel and, from time to time, to swap the current microstate 
(i.e. the instantaneous set of atomic positions and momenta) 
between two parallel runs. The swapping is accepted/rejected 
according to a Metropolis Monte Carlo criterion, namely with 
the probability 

p =min { 1 . (5) 

exp[(V V ({x, p} ft ) - V V ({x, p} p .)) (Pt - Pj) } 

where V^({x,p}jg.) and V^^({x,p}^.) are the potential en- 
ergies of the two microstates, respectively at ft = £#7] and 
/3 j = ksTj in the phase space points {x,p}jg. and {x,p}^. at 
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the moment of the attempted swapping. If the swap is ac- 
cepted, then the microstates corresponding to temperatures 7] 
and Tj are simply interchanged. If the swap is rejected, the 
microstates are further aged at their own temperature. The 
key feature of this method is that the sampling of the system 
phase space obtained by the piece-like replica exchange tra- 
jectories is consistent with the canonical (conditional) proba- 
bility density function at each target temperature. However, 
since the individual pieces of the replica exchange trajecto- 
ries are obtained by swapping from higher temperatures, they 
more likely overcome possible free energy barriers. In short: 
the replica exchange trajectories are ergodic. 

The simulation techniques described above (and the collec- 
tive variables presented in the next section) have been im- 
plemented in the CMPTool simulation package! 18 n 2Q l i n par- 
ticular, the combination of the restrained MD with Parallel 
Replica technique allows a two-level parallel approach. The 
fist level of parallelism is over the replicas while the second, 
implemented according to the usual domain decomposition 
paradigm, is within each replica. While the latter level of par- 
allelism, that requires a tight connection, was implemented 
using the Message Passing Interface (MPI) application pro- 
gramming interface, the former was implemented at a script- 
ing level. This allowed to run the simulations on a cluster 
of loosely connected multicore machines communicating by a 
Gigabit network. 



above definition of the collective variable by a smooth analyt- 
ical approximation of it that, in a proper limit, converges to its 
exact definition and performed biased MD runs according to 
this representation of 3%{s). This smooth analytical approxi- 
mation is obtained in two steps: (i) first we obtain an analytic 
(explicit) expression of mhi/{|x c — xf\} as a function of the 
positions xf, and (ii) then we introduce a smooth approxima- 
tion to this expression. The first step consists in recognizing 
the following identity: 
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in{|x c -x?|} = £ 
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(6) 



where G (x) is the Heaviside step function. Let / be the O atom 
closest to the center of the nano-particle, then n^ ®(l x c — 
xp| — |x c — x?|) = Su , where 5// is the Kronecker symbol. 

This implies that the result of this sum is |x c — xf\, i.e. the 
distance from x c of the closest O atom. Eq. [6] is, therefore, 
the analytical expression of the collective coordinate «^(x). A 
smooth approximation to &(x) can be obtained by replacing 
the Heaviside step function by a sigmoid function. We used a 
sigmoid function expressed in term of the Fermi function: 



B. Collective variables for modeling the crystallization process 
in nanoparticles 

We now discuss the collective variables used to study the 
crystallization process. It is worth to stress again that in this 
paper we are not interested in studying the detailed mecha- 
nism of crystallization in nanoparticles. Rather, we investi- 
gate the "relative stability" of the ordered vs the disordered 
phase as a function of the size of the nanoparicles and the tem- 
perature. Therefore, our collective variables must be able to 
distinguish between the crystalline and the amorphous phase 
(i.e. they need to be order parameters), rather than modeling 
the mechanism of the crystallization. 



1. &(x) order parameter to monitor the size of the nanoparticle 

We introduce the notion of size of the nanoparticle by a col- 
lective variable denoted by the symbol &(x). In a system in 
which the nanoparticle is made of atoms of type 'A' (Si in this 
case) and the matrix is made, or contains, atoms of type 'B' (O 
in this case) a possible definition of the collective coordinate 
&(x) is the distance between the center x c of the nano-particle 
(a point kept fixed during the simulations) and the closest 
oxygen atom, i.e. &(x) = mim{|x c — xf\} 9 where xf is the 
coordinate of the i-th oxygen atom. The force acting on the 
atoms associated to this collective variable cannot be straight- 
forwardly evaluated since M(x) is a non-analytical function 
of x and, therefore, there is no way to proceed through the 
direct calculation of Vk f /2(&(x) -&*) 2 . We replaced the 
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where A is the parameter controlling its smoothness. In our 
simulations A has been chosen such that the sigmoid function 
goes from 0.95 to 0.05 in one atomic layer (~ 0.2 nm). A 
consequence of this replacement is that the size of the nano- 
particle is now defined as a weighted average of the distance 
of one atomic layer of oxygen atoms from the centre of the 
nano-particle. 



2. ( x ) order parameter to monitor the phase of the 
nanoparticle 

We compute the free energy of a Si nano-particle embedded 
in silica as a function of its degree of order, as measured by 
the bond-orientational order parameter (^(x)) introduced by 
Steinhardt et al. 21 for bulk systems. In this paper we adapted 
the original definition to the case of confined systems, as de- 
scribed below in detail. 

In general, £?e(x) is defined as 



M*) = [j— - £ |^ 6m (x)| 2 (8) 

\ * ra=— 6 / 

where i?6m ( x ) is the normalized and weighted sum of atomic 
vectors q l 6m (x) (defined below) which, in bulk systems, reads 
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^6m(x) 



(9) 



where N is the number of atoms in the system, Ni is the num- 
ber of nearest neighbors of the atom i and m = —6, .... 6. In 
the case of confined systems we limit the sum over i to just 
the atoms belonging to the nano-particle. Consistently with 
our definition of the size of the nanoparticle, we identify these 
atoms as those at a distance lower than ^* from the center 
of th e nano particle (0$* is the size of the nanoparticle, see 
II B l\ . According to this definition, the ^m(x) of the 
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nanoparticle is: 



£f =1 ^ m W(l-0(|xf-^|-^*)) 
Hi* 



(10) 



where 0(x) is the Heaviside step function. As for the case 
of the collective variable ^(x), the Heaviside step function 
is, in practice, replaced by a sigmoid function S(t) = 1 — 
l/(l+exp[Ar]). 



The q l 6m (x) function appearing in Eq. 
ing to the following expression: 
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^6m( x ) 



is defined accord- 



(11) 



where Y^m^ij) is the spherical harmonic function of degree 6 
and the component m computed on the solid angle x\j formed 
by the distance vector x\j and the reference system. The sum 
runs over the Ni nearest neighbors of the atom i. 

The sum over the m component in Eq. [5] makes the collec- 
tive coordinate ^(x) rotationally invariant, i.e. independent 
on the orientation of the reference system. 

When the system is an ideal crystal and the temperature is 
K, the environment of all the atoms is the same and, there- 
fore, i?6( x ) is maximum as there is not interference among 
the q l 6m (x). On the contrary, in a perfectly disordered sys- 
tem the orientation of bonds is random and, therefore, there is 
complete interference among the q l 6m (x), and ^(x) is zero. 
However, even when the system is at finite temperature and 
its size is finite, this order parameter is still able to distinguish 
between a disordered and a crystalline phase. 

Before concluding this section it is worth to mention 
that the ^ collective variable, or its modifications, has 
been already used to study crystallization by atomistic 
simulations^ and experiments P3 



C. Simulation setup 

In the present investigation, the restrained MD is governed 
by the superposition of the physical potential and the restrain- 
ing potential k/2(£ 6 (x) -£*) 2 -f £'/2(^(x) -^*) 2 . kandk! 
are the parameters controlling the degree of biasing and must 
be large enough such that along the MD the values of £?6 (x) 



and &(x) oscillate about the target values and respec- 
tively. In this work we use the Billeter et alP3 environment- 
dependent force field. This classical force field, which is an 
extension of the Tersoff potential, is defined as the sum of 
generalized Morse pair potential terms: V(x) = Y,i>j v u(\ x i ~ 
Xj\), where / and / denote the chemical species of atoms 
i and j, respectively. The pair potential vjj{x) is given by 
v/jW =///M[A//exp[-A//*] -b I j(x)A I jQxp[-fi I jx}}. f u (x) 
is a switching function that makes vjj(x) to go smoothly to 
zero at the cutoff distance rfj 1 (rfj* is such that v/j(jc) is non 
zero only between nearest neighbor atoms). bu(x) is the en- 
vironment dependent term, which is function of an effective 
coordination number. The effective coordination number em- 
bodies three-body terms through the angle formed by i, j and 
all their nearest neighbors. In addition to the terms mentioned 
above, the Billeter et al. force field contain terms that pre- 
vent unphysical over/under-coordination at the Si/SiC>2 inter- 
face. For a detailed description of the potential we refer the 
reader to the original paper. The reliability of this force field 
in modeling equilibrium and dynamical properties of Si nano- 
particles embedded in silica, of the Si/a -SiO? interface and Si 
nanowires has been already established! 2 6 l 28 n 3 1 1 

The computational samples are prepared by thermally an- 
nealing a periodically -repeated amorphous silica system and 
embedding Si nano-grains (extracted from a well equilibrates 
either amorphous or crystalline bulk). The amorphous sil- 
ica sample is prepared through the quenching-from-the-melt 
procedure, that is by cooling down very slowly a high tem- 
perature SiC>2 melt. The total system contains from ~ 6000 
to ~ 12000 particles, corresponding to a nano-particle radius 
varying in the range 1—2 nm. Computational samples are 
first thermalized at 300 K and ambient pressure by using the 
Martyna-Tobias-Klein variable cell algorithm®! in order to re- 
lease possible stress at the Si/silica interface. Typically, during 
such a thermalization step, the nano-particles slightly shrink. 
After this initial step, we impose the restraint on the size of 
the nano-particles and thermalize the samples at the various 
target temperatures. Because of the restraint on their size, at 
this stage the size of the nano-particles does not change. After 
this treatment the samples are ready for the restrained simu- 
lations described above. The simulations are performed at a 
fixed volume. However, we checked that the pressure is close 
to the ambient value all along the simulation. 

In order to verify possible artifacts due to finite-size effects, 
we repeated the calculation of the observable < k(0(x) — 
0*) >H kk , at few selected value of J?^ and ^* on samples 
of different size of the silica matrix. We did not observe any 
significant difference in them (the differences were within the 
statistical error). This demonstrates that there are no finite- 
size effects in our free energy calculations. 

We compute the free energy F(£?6,&) in the range JS^ G 
[0,0.65] and Si e [0.8, 1.8] nm. The rationale for the upper 
limit of the range is that the value of the bond-orientational 
order parameter for an ideal Si crystal is J?6 = 0-63 and, since 
from experiments and previous calculations it is known that 
Si crystalline nanoparticles assume a structure with a (dis- 
torted) diamond-like core and a disordered peripher^ 1 ^ 7 ^ 8 ^ 33 ! , 
we expect the i?6 of crystalline nanoparticles be lower than 
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this limit. The samples created according to the protocol de- 
scribed above confirm that the JS^ of crystalline nanoparti- 
cles is lower than this upper limit. However, after a prelim- 
inary scanning of V^ 6 F(^^) that allowed to identify the 
region of JS^ containing the minima in the above domain, we 
restricted the calculations to a smaller range: [0.04,0.285], 
[0.07,0.35] and [0.03,0.38] for the nanoparticles of radius 0.8, 
1.3 and 1.8 nm, respectively. As for the size of nanoparticles, 
the same experiments mentioned above report unusual phase 
transitions (disorder-to-order with growing T) for nanoparti- 
cles in the range [0.5,2.0] nm. We decided therefore to study 
nanoparticles of three size in this range: 0.8, 1.3 and 1.8 nm. 

Finally, the restrained MD/parallel tempering simulations 
where performed at 300, 500, 750, 1000, 1250 1500 and 
1750 K. 



III. RESULTS AND DISCUSSION 

As a first preliminary step, we begin the presentation and 
discussion of our results by motivating the use of the restraint 
method in combination with the replica exchange method 
through a relevant example. In Fig. [T] two quasi-crystalline 
configurations are shown, corresponding to the same value of 
the parameter. These configurations embed different de- 
fected structures. The configuration shown in the left panel is 
characterized by an extended disordered region in the bottom- 
right part of the nano-particle. At a variance, two smaller dis- 
ordered regions characterize the second configuration shown 
in right panel, respectively in the bottom-right and top-left part 
of the nano-particle. Both configurations should be consid- 
ered for the correct evaluation of the gradient of the free en- 
ergy < k(6(x) — 0*) >H kh! corresponding to the same value 
of £?l . However, a sizable free energy barrier likely separates 
them. Therefore, if the simulation is started from one of the 
two configurations, then the second one likely could not be 
visited during the explored time scale. By using the replica 
exchange method we are able to properly and efficiently sam- 
ple both of them and compute dF(0)/d0 at the present 
accurately. 

As a second preliminary step, we summarize the experi- 
mental findings we aim to address. By comparing Energy Fil- 
tered Transmission Electron Microscopy (EFTEM) and Dark- 
Field Transmission Electron Microscopy (DFTEM) images in 
Si-rich SiO x samples it was shown that Si nano-particles start 
to form at 1000°C. 10 At this temperature they are all amor- 
phous, while at 1100°C about one third become crystalline. 
By further increasing the annealing temperature by 50° C, the 
fraction of crystalline nano-particles rises up to 60%, while 
the average size of the nano-particles and the distribution of 
their size remains almost unchanged. Finally, at the annealing 
temperature of 1250°C 100% of nano-particles are crystalline. 
At this temperature the average size is slightly increased, but 
the particle size distribution is still largely superimposed to 
the distributions observed at 1100°C and 1150°C. It was also 
found that the system is stationary with respect to the amor- 
phous vs. crystalline composition. Similar investigations have 
been performed on Si/SiC>2 multilayers 10 where the growth of 




Figure 1 . Two different configurations of an embedded silicon nano- 
particle with radius as large as 1.8 nm. They both correspond to 
=0.19. Oxygen atoms are displayed in red and Silicon atoms 
in yellow. In order to improve the readability, only the atoms laying 
within a 1.5 nm- thick slice are drawn. 



the crystalline fraction with the annealing temperature is even 
more sudden: the degree of crystallinity increases from about 
15% to 90% when the annealing temperature is increased from 
1100°C to 1200°C. Also in this case it was demonstrated that 
the samples are at the equilibrium. 

We can now turn to the results of our simulations. In 
Fig. [2] the free energy curves F{£1§) of Si nano-particles of 
size^* = 0.8 nm, ^ = 13 nm, and = 1.8 nm at var- 
ious temperatures in the range 227° C - 1477°C are shown, 
as obtained from our simulations (our calculations are per- 
formed in Kelvin, while the results are presented in Celsius 
for homogeneity with available experimental data). Above 
and throughout the paper we shall denote the radius of the 
nanoparticle by the symbol which is the target value of 
the collective variable and, via the corresponding re- 

strain term (see Sec. |IIQ , defines its size. As a first remark, 
we notice that two metastable states, one a low J^6 and one 
at higher £1$, are present at all temperatures and sizes. In 
general, to high values of correspond crystalline states 
while to low values of the same order parameter correspond 
disordered (amorphous) states. However, especially for the 
smallest nanoparticle, the difference between the value of J?6 
corresponding to the two metastable states is small. There- 
fore, the identification of the phase corresponding to a given 
state requires a further investigation of the structure. We per- 
formed this analysis by computing the Si-Si partial pair corre- 
lation function g(R) for the atoms belonging to the nanopar- 
ticle (|xf — x c | < ^*) on the configurations corresponding 
to the two metastable states (£?e(x) = where ^ corre- 
sponds to one of the two minima). In the left-most panel of 
Fig. [3] we report the g(R) of the low (bottom panel) and high 
(top panel) metastable states for the 8% = 0.8 nm nanoparticle 
at various temperatures. For the sake of comparison, we also 
show the g(R) of the bulk amorphous and crystalline states at 
T = 627°C (i.e. T = 900K). For the high £! 6 metastable state, 
we notice that even at the highest T the g(R) is characterized 
by three peaks in the range < r < 5. These peaks correspond 
to the bulk-like first, second, and third neighboring shell, re- 
spectivelly. They broaden and reduce in height by increasing 
the temperature; nevertheless, they are still well visible even 
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Figure 2. Free energy vs J2(, curves for nano-particles with radius 0.8 nm, 1.3 nm and 1.8 nm. The curves are shifted to improve readability. 



at T — 1227° C. In general, even at low T, these peaks are 
broader than the corresponding bulk crystalline ones. This is 
due to a structure of the nanoparticle which is composed of 
two regions : a crystal-like core and a disordered surround- 
ing shell. The latter is in contact with the S1O2 matrix. This 
two-region structure is consistent with the structure found by 
Hadjisavvas and Kelires 33 in their investigation on crystalline 
nanoparticles embedded in silica. We shall provide further 
evidence of this two-region structure below. 

Let us move to the analysis of the structure of the low J?6 
metastable state. As a first remark we notice than in the same 
R range analyzed for the high £?6 case there are only two 
peaks. The first one, sharp and intense, corresponds to the 
nearest neighbor Si-Si pairs. The second one, very broad and 
small, is usually assigned to the second neighbor pairs, which 
in amorphous system are distributed over a large r range. 
There is no other peak in the < r < 5 domain. Compar- 
ing the g(R) of the low metastable state of this particle 
with amorphous bulk Si we notice that there is a one to one 
correspondence between the equivalent peaks in the two sys- 
tem. Similar results, both for the low and high J?6 metastable 
states, are found for the £& — 1 .3 nm and ^=1.8 nm nanopar- 
ticles (see the central and right-most panels of Figs. [3j respec- 
tively). On the basis of this analysis, we identified the high J?6 
metastable state of all the nanoparticles at all temperatures to 
be of crystalline nature while the one at low to be of amor- 
phous type. 

The above conclusions are confirmed by the visual inspec- 
tion of the structure of the nanoparticles in the low and high 
^6 domains. In Fig.|4]we show few snapshots collected along 
the restrained MD at the values of the £?6 collective variable 
corresponding to the minimum of the free energy in the amor- 
phous and crystalline metastable states, respectively, at low 
and high temperatures. It appears evident that the nanoparti- 
cles in the high £?6 domain have a crystal-like core in which 
the tetrahedral orientation of the Si- Si bonds is preserved. 
This core is surrounded by a shell containing disordered re- 
gions. At a variance from this, the structure of the nanopar- 
ticles corresponding to the low £?6 free energy minimum are 
completely disordered. 

We further investigated the structure of the crystalline 
metastable states by computing the shell-by-shell J?6- This 
quantity, denoted by the symbol £}(>(R), is obtained by lim- 
iting the sums in Eq. [10| to the atoms laying in the spherical 



layer of thickness 0.1 nm at the distance R from the centre 
of the nanoparticle. In Fig. [5] is shown the <S^(R) for the 
three nanoparticles at the same T as in Fig. [4] As a first re- 
mark, Fig. [5] confirms the observation that the degree of or- 
der decreases in going from the center to the periphery of 
the nanoparticle. This is consistent with the results of pre- 
vious works .El However, Fig. [5] also indicates that there is a 
qualitative difference in the structure of nanoparticles of dif- 
ferent size. For the nanoparticle of size 0.8 nm we notice 
that the J2e(R) decreases monotonically with R. The trend 
is very similar both at low and high T. On the contrary, al- 
ready for the nanoparticle of size 1.3 nm, we observe that the 
J£(){R) is characterized by one plateau region in the domain 
nm < R < 0.5 nm. Beyond this region, the i?6 quickly de- 
creases, till reaching the value of ~ 0.3 at the Si/a-Si02 inter- 
face. For the largest nanoparticles, at low T we observe two 
plateau regions, one in the domain nm < R < 0.5 nm and 
another, at lower in the domain 0.5 nm < R < 1.1 nm. 
At higher T we have only one plateau region in the domain 
nm< R < 1.1 nm. In both cases, beyond R = 1.1 nm the 
^6 g° es ver Y quickly to the value of ~ 0.35. These observa- 
tions indicate that there is a qualitative difference between the 
structure of the crystalline nanoparticle with their size, namely 
that there is a threshold below which the ordered nanoparticle 
does not have a crystalline core with a homogeneous degree 
of order. This difference is preserved also at higher T 

Once having identified the nature of the low and high J?6 
metastable states, and having analyzed their structure, we turn 
to the analysis of the order-disorder phase transition with the 
size of the nanoparticle and the temperature of the system. 
Our simulations (see Fig. [2]) provide the following qualitative 
sharp picture: for small nano-particles (^* =0.8 — 1.3 nm) 
at low temperature (T < 727° C) the most stable configura- 
tion corresponds to the amorphous phase, while the crystalline 
state is found to be more stable at higher temperatures. On 
the contrary, for larger particles (ffl* > 1.8 nm) this behavior 
is inverted, resulting similar to bulk- like conditions: at low 
temperatures (T < 977° C) the crystalline phase is the most 
stable one, while the amorphous phase is preferred at higher 
temperatures. Interestingly enough, for small nano-particles 
the equilibrium temperature (i.e. the temperature at which the 
free energy of the disordered and ordered phase are the same) 
decreases with the increase of the size of the nano-particle. 
This is indeed an effect of the steady increase of stability of 
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Figure 3. Si-Si pair correlation function (g(R)) of the low and high J2(, metastable states at various T. Bulk crystalline and amorphous g(R) 
are also reported for comparison. 



the crystalline phase with respect to the disordered one with 
the size of the nano-particle. The results of our simulations 
provide the following picture, consistent with the experimen- 
tal results: 5 10 at low annealing temperature the nano-particles 
are small and amorphous as, due to the inversion of stability 
with respect to bulk-like systems, this is thermodynamically 
the most stable phase; at moderately higher temperatures the 
average size and the size distribution of the nano-particles are 
slightly changed (the average size is slightly increased). The 
largest particles in the sample transform from amorphous to 
crystalline, the most stable phase for large nanoparticles at this 
T. By further increasing the temperature the average size of 
the nano-particles increases significantly and the larger ones 
tend toward the crystalline state (i.e. they follow the change 
in stability from disorder to order, as induced by their grow- 
ing size). On the other hand, the smaller particles undergo a 
amorphous-to-crystalline transition due to the increase of the 
temperature and the inversion of the stability with respect to 
the bulk-like system. Even in this case they eventually crys- 
tallize. The model above is based on the experimental obser- 
vation of the dependency of the average size and size distri- 
bution of the nanoparticles on the temperature. 10 Our results 
bring to the conclusion that the observed crystallization of the 
nanoparticles with the increase of the temperature is due to the 
interplay of the effect of the temperature on their size and the 
inversion of the relative stability of the amorphous and crys- 
talline phase with the temperature for small nanoparticles. 

IV. CONSLUSIONS 

In this paper we investigated the relative stability of the 
amorphous vs crystalline nanoparticles of size ranging be- 



tween 0.8 and 1.8 nm. We found that, at variance from bulk 
systems, at low T small nanoparticles are amorphous and 
they undergo to an amorphous-to-crystalline phase transition 
at high T. On the contrary, large nanoparticles recover the 
bulk-like behavior: crystalline at low T and amorphous at high 
T. 

We also investigated the structure of the crystalline 
nanoparticle. Our results, in agreement with previous 
works, 33 demonstrate that this kind of nanoparticle are formed 
by an ordered core surrounded by a disordered periphery. 
However, they also indicate that the details of the structure 
of the crystalline core depend on the size of the nanoparticle 
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